POWER LAW STATISTICS OF AVALANCHES IN MARTENSITIC 

TRANSFORMATION 
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We devise a two dimensional model that mimics the recently observed power law distributions 
for the amplitudes and durations of the acoustic emission signals observed during martensitic trans- 
formation [ Vives et al, Phys. Rev. Lett. 72, 1694 (1994)]. We include a threshold mechanism 
arising from the athermal nature of transformation, long-range interaction between the transformed 
domains, inertial effects, and dissipation arising due to the motion of the interface. The model 
exhibits thermal hysteresis of the transformation, and more importantly, it shows that the energy is 
released in the form of avalanches with power law distributions for their amplitudes and durations. 
Computer simulations also reveal morphological features similar to those observed in real systems. 
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Many spatially extended driven systems naturally 
evolve to a marginally stable state characterized by 
avalanches with power law distributions for their am- 
plitudes and durations reflecting lack of intrinsic length 
scales and time scales in the system. Such a state is 
termed as self-organized critical (SOC) state by Bak et 
al 0| . Several physical systems exhibit SOC features; for 
example, earthquakes ||| , acoustic emission from volcanic 
rocks ||, and stress drops during jerky flow Q, to name 
a few. Recently, Vives et al ||, measured the acoustic 
emission (AE) signals during martensitic transformation 
of Cu-Zn-Al single crystals under thermal cycling. They 
reported power law statistics for the amplitudes and du- 
rations of the AE signals both during cooling and heating 
runs. To the best of our knowledge, there is no strain (or 
displacement) based model of martensitic transformation 
which explains these results. More over, any prospec- 
tive model has to take into account the non-equilibrium 
nature of the hysteresis. Even though extensive theo- 
retical studies exist on martensitic transformations |i] ^) , 
the influence of dissipation and defects on hysteresis has 
received very little attention. Here, we propose a sim- 
ple phenomenological model which captures the power 
law distribution of AE signals along with the thermal 
hysteresis of the transformation. Below, we will briefly 
collect SOC type dynamical features of the martensitic 
transformation relevant for modelling the system. 

Martensitic transformation is a first-order, solid-solid, 
diffusionless, structural phase transition. On cooling, the 
unit cell gets distorted leading to the nucleation of 
thin plate- like product domains with a twinned structure. 
This induces internal strains, which in turn induce long- 
range fields, that block the transformation leaving the 
system in a two phase metastable state. ( Note that this 
implies the existence of built-in threshold mechanism.) 
Thus, the amount of the transformed phase is entirely 
determined by the excess free energy and, an additional 
undercooling is required for further growth. This implies 
that thermal fluctuations have little role in the transfor- 
mation kinetics. Thus, the transformation is athermal 
and hence the nucleation is athermal, usually occurring 
at defects like dislocations P,H-fO|- Further, the emis- 



sion of acoustic energy in the form of bursts implies that 
inertial effects are important. Indeed, Bales and Good- 
ing, and later Reid and Gooding Q, have studied the 
importance of including the inertial effects. Vives et al 
attribute the mechanism of irreversible release of the elas- 
tic energy in the form of avalanches, each of which cor- 
respond to the motion of one (or more) interface, to the 
evolution of the system from one metastable state to an- 
other |). Since the interface moves at near velocity of 
sound, there is dissipation associated with it. We include 
all these features along with the long-range interaction 
between the transformed domains. 

For the sake of simplicity, we consider, a 2d square 
to rectangle transition, for which the free energy is usu- 
ally a function of all the three components of strain de- 
fined by ei = (r) xx + r) yy )/V2, e 2 = (j] xx - % y )/V2, and 
e3 = r\x y = Vyx, where ?/y = |(§^ + g^.) refers to the 
components of the strain tensor and m is the displace- 
ment field in the direction i ( i = x,y). Here, e\ is the 
bulk dilatational strain, e 2 is the deviatoric strain and e% 
is the shear strain. In athermal martensites, deviatoric 
shear strains play a dominant role in the transformation 
kinetics Ifj). In this paper, we consider the deviatoric 
strain e 2 = e(r) as the principal order parameter since 
volume changes are small The effect of other 

components of the strain, i.e., the bulk and shear strain 
is accounted phenomenologically by considering a long- 
range interaction between the deviatoric strains. 

We write the rescaled free-energy functional as 
F{e{r)} = F L {e{r)} + F lr {e{f)}, where F L is 



F L = / df 



(1) 



where D and a are in a scaled form. Fi r is an effec- 
tive long-range term that describes transformation in- 
duced strain-strain interactions. In Eqn. (1), fi(e(r)) = 
^e(f) 2 — e(r) 4 + is the usual Landau polynomial 

for a first order transition, where ST = (T-T C )/(T -T C ) 
is the scaled temperature. T is the first-order transition 
temperature at which the free energy for the product and 
parent phase are equal, and T c is the temperature be- 
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low which there are only two degenerate global minima 
e = ±6m- Following Cao et al JL(|, the effect of localized 
defects acting as nucleation sites is simulated by an in- 
homogeneous stress field, o~(r). This term modifies the 
free-energy // in a way that renders the austenitic phase 
locally unstable leading to the nucleation of the prod- 
uct phase. As mentioned in the introduction, the phys- 
ical cause of the long-range interaction is the coherency 
strain between the parent and the product phase. Such 
an interaction is also expected to arise due to the cou- 
pling of e(r) with the other components of the strain 
order parameter. Recently, an effective long-range inter- 
action has been shown to result, both for the bulk [ p~3|Jl4| 
and interface between the deviatoric strains if e[ 
and e§, are eliminated by imposing (St. Venant) elas- 
tic compatibility constraint. In ref H], to describe the 
effect of the interface between the austenitic-martensite 
phases, the authors introduce an explicit interface term. 
Instead, we follow Wang and Khachaturyan [|l5| who have 
shown that the interface can be naturally described by 
accounting for coherency strains at the parent-product 
interface by a symmetry allowed fourth order anisotropic 
long-range interaction in the free energy. Following this, 
we write down the long-range term phenomenologically 
in the Fourier space as 



F lr {e}= / dkB(k/k){e 2 (r)} k {e 2 (f)} k „ 



(2) 



where {e 2 (j)}fc and {£ 2 {r)} k » are the Fourier transform 
of e 2 (r) and its complex conjugate respectively. For 
the square to rectangle transformation, the favourable 
growth directions for the product phase are the habit 
plane directions [11] and [11]. In addition, free-energy 
barriers should be large along the [10] and [01] direc- 
tions. These features are well captured by the simple 
kernel B(k/k) — -~^f39(k — A)k 2 k 2 , where k x and k y are 
the unit vectors in the x and y directions, and /3 is the 
strength of interaction and A in the step function 9(k—A) 
is a cutoff on the interaction range. This kernel incorpo- 
rates the effect of the interface in a natural way since 
the cost of growth progressively increases with the trans- 
formation in directions where the kernel is positive which 
not only aids growth along the habit plane directions but 
also limits the growth of domains transverse to it. We 
stress that this is only a simple choice and is not unique. 
Other kernels with similar orientation dependence will 
give similar results. The real space picture of B(k/k) is 
similar to the long-range interaction of Kartha et al [[13| . 

Even though deviatoric strain is the order parameter, 
basic variables are the displacement fields. Thus, we start 
with the Lagrangian L = T — JF, where F is the total 
free-energy with the kinetic energy T defined by 



T 



drp 



du x (r,t) 
dt 



du v (r,t) 
dt 



(3) 



Here p is the mass density. The dissipation associated 
with the movement of the parent-product interface is 



represented by the Rayleigh dissipative functional p6[ . 
Further, since deviatoric strains are the dominant ones, 
dissipative functional is written entirely in terms of e(r), 

ie., R — \l j dr(-gte(r,t)^ . (This is consistent with the 
fact that shear and bulk strains are known to equilibrate 
rapidly and hence do not contribute.) The equations of 
motion for the displacement fields are 
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Using the above equations, we obtain the equation of mo- 
tion for e(r, t), which after scaling out p and D (in terms 
of rescaled space and time variables) gives 



dt 2 



e(f,t) 



SF d ,„ ; 



(5) 



Thus, both (3 and 7, here, are to be taken as rescaled 
parameters. The structure of Eqn. (5) is similar to that 
derived in |ll|] for 1-d except for the long-range term. 

Simulations are carried out by discretizing Eq.(5) on a 
Nx N grid using the Euler's scheme with periodic bound- 
ary conditions. The mesh size of the grid is Ax = 1 and 
the time step At — 0.002. The long-range term is com- 
puted using psuedospectral method: this term is calcu- 
lated in the Fourier space and then inverse fast Fourier 
transform is effected to obtain the real space interaction. 
(We have used the cutoff distance A = 0.2.) We have 
studied both the single defect and many defect nucleation 
case. Here, we reports results on multi-site case and only 
mention results on single site case wherever necessary. 

Consider the nucleation and growth for a single 
quench. We use a random distribution of defects and rep- 
resent their stress field by er(r) = 53j ax a o{^j) ex P[~ 
r j| 2 /C 2 ]i where r*j refers to the scaled coordinates of the 
defect sites and j m ax is the total number of defect sites. 
uo(fj) is taken to be uniformly distributed in the inter- 
val [—0.3,0.3]. Initially, the system is taken to be in 
a homogeneous state with e(r, 0) distributed uniformly 
in the interval [-0.005,0.005]. At t = 0, we switch on 
the stress-field a(r). Figure 1 shows snap shots of the 
transformation at t = 10, 12, 15 and 50. (Grey regions 
represent the austenitic phase e = 0, black and white 
regions represent the two variants.) Nucleation of the 
product phase occurs with a value e ~ ±2 at several 
locations where the magnitude of the stress-field is suffi- 
cient to make the system locally unstable ( t = 10). We 
note here that the value of e for the two variants is larger 
than that given by fi(e(r)) alone ( cm ~ ±1.31 for these 
parameters) due to the long-range interaction. In a short 
time, we see the other variant being created adjacent to 
these nuclei. By t = 12 ( Fig. 1), the structure further 
develops into twinned arrays, propagating along [11] and 
[11] directions respectively. We also see that several new 
domains are nucleated at a finite distance from the preex- 
isting domains which is consistent with the autocatalytic 
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nucleation mechanism known to operate in martensitic 
transformation J|||. Though, the new nucleation sites, 
most often coincide with the defect sites, occasionally it 
is seen in defect free regions. ( Similar observation was 
made in the single site case also.) This can be attributed 
to the stress accumulation at these sites as a result of 
the long-range term due to the preexisting martensitic 
domains. Note that the twinning is irregular which is 
again due to the mutual interaction between the vari- 
ous domains. There is very little growth beyond t = 30. 
Thin needle-like structures can also be seen to emerge 
from larger domains ( t — 50) as reported by Vives et al. 
There is also a distribution of domains sizes. 

We simulate thermal cycling of the transformation by 
continuously changing ST from + 40 to -80, and back in 
a duration of t = 1000 units at a constant rate, both for 
the heating and cooling runs. For the reverse cycle, the 
final configuration of the cooling run is used as the ini- 
tial configuration. Figure 2 shows the area fraction <pA 
versus ST for the cooling and heating runs ( o ). In the 
cooling run, transformation starts around ST rns ~ —2.0 
showing a rapid increase in (f>A up to ~ 30%, thereafter, 
it increases linearly up to 90%. The transformation is 
complete by ST m f ~ —59. In the heating run, the parent 
phase appears only at ST as ~ —26.0 and 4>a decreases 
almost linearly till the transformation is nearly complete 
around ST a f ~ 18. ( For the sake of comparison, we have 
also shown the single-site hysteresis loop by •. ) 

The most important feature of the model is that the 
changes in <pA are actually jerky which can only be seen 
on a finer scale. Since, the rate of energy dissipated 
dE/dt = -2R(t), we have calculated R(t) ( or R(ST)). 
In Fig. 3, we have plotted R(t) with the inset show- 
ing the enlarged section of the peak which clearly shows 
that the rate of energy release occurs in bursts consistent 



(Latent heat also 
Further, the dis- 



with the acoustic emission studies 
shows a pattern similar to Fig. 3 
tribution of the amplitudes of R(t) denoted by Ra has a 
tendency to approach a power law Dr(Ra) ~ i?^ QH with 
an ~ 2.6 (o in Fig. 4. We have also shown the single-site 
results by •.) Similarly, the distribution D R (At) of the 
durations At of the energy bursts scales with At given 
by D R {At) ~ At~ TR with tr ~ 3.2. Although, the scal- 
ing regime is almost identical to Dr(Ra), we find that 
the scatter is slightly more in this case. ( This is true 
of experimental results as well|5| and in SOC models 
0). The conditional average |l8] of Ra, for a given At 
denoted by < Ra > c behaves as < Ra) > c ~ At XR . We 
get the exponent value xr ~ 1.36. We also find that the 
scaling relation tr = xr(oir — 1) + 1 is satisfied quite 
well. In experiments, however, one measures the am- 
plitude of the AE signals A ae , ie., Ra ~ A^ e . Using 
the relationship between the two joint probability distri- 
butions D(RA,At) oc D(A ae , At)/A ae , it can be easily 
shown that ctR = (a ae + l)/2 where a ae is the expo- 
nent corresponding to A ae . The other two exponents 
remain unchanged. Using the experimental values || ( 
a ae ~ 3.8, r ae ~ 3.6 and x ae ~ 1), we see that or ~ 2.4. 



Thus, we see that our values are in reasonable agreement 
with experiments considering the fact that real systems 
are 3-d. It must be stated here that in 3-d even the num- 
ber of variants are generally more and one also expects 
that the mechanisms operating in 3-d cannot be fully ac- 
counted for in 2-d |8p5[| . We have also carried out a sim- 
ilar analysis on R(t) for the heating run. Even though, 
the changes in R(t) occurs in bursts in the heating run as 
well, we find that the scatter in the distributions is more 
than that for the cooling run. 

In conclusion, the fact that the energy release occurs 
in the form of bursts with a power law statistics for 
the avalanches is well captured by the model. In addi- 
tion, the model shows hysteresis under thermal cycling. 
The power law statistics can be attributed to the fact 
that we have included important ingredients of SOC dy- 
namics, namely, the threshold dynamics, dissipation, the 
generation of large number of metastable states during 
the transformation, and a relaxation mechanism for the 
stored energy. This relaxation time scale corresponds to 
the inertial time scale as can be inferred from the fact 
that the basic variables are the displacement fields which 
in turn set the limit on the fastest time scale. Indeed, 
from our simulations we find that the interface move- 
ment occurs at time scales of a few units of (scaled) time. 
Compared to this the driving force increases with tem- 
perature at a slow pace. The important feature of in- 
ducing large number metastable states during cooling or 
heating runs is due to the long-range interaction between 
the transformed domains as can be seen from the follow- 
ing reasoning. We note that the value of this term at any 
location is the result of the superposition of the contribu- 
tions arising from the spatial distribution of the already 
transformed domains which in turn leads to a complex 
terrain of local barriers ( metastable states). These self 
generated ( transformation induced) local thresholds, at 
a given time, must be overcome by the increase in the 
driving force arising from the slow cooling (or heating). 
We note that once a local barrier is overcome, part of 
the driving force goes in creating a new twin and the 
rest is dissipated in the form of burst of energy ( due 
to the advancing interface) . The fact that long-range in- 
teraction is at the root of creating the local thresholds 
is further supported by the fact that we find a power law 
distributions even in the single site nucleation case. (See 
• in Fig. 4.) The presence of defect sites only triggers 
the initial nucleation process. This must be contrasted 
with disorder based Ising models |HJ which also produce 
power law statistics for avalaches and field induced hys- 
teresis. However, by subjecting the samples to repeated 
thermal cycling, Vives et al have verified that in marten- 
site transformation, it is the dynamical ( transformation 
induced) disorder that is at the root of the avalanches. In 
this sense, our model is the first to capture both thermal 
hysteresis and jerky nature of the transformation based 
on dynamical disorder and is independent of quenched 
disorder. Finally, the morphological features, including 
the needle shaped domains are very similar to those ob- 
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served in experiments. 
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Figure captions: 

Fig. 1: Snapshots of the morphology of the transfor- 
mation at t = 10, 12, 15 and 50 time steps for multi-site 
nucleation. The parameters values are: N = 128, ST — 
-2.0, (3 = 50, 7 = 4, ct - [-0.3, 0.3] and j max = 164. 

Fig. 2: Area fraction (f>A as a function t for cooling 
and heating (o) runs for N = 256, (3 = 50,7 = 4, <j = 



[-0.3,0.3] and jmax ~ 1% of iV 2 . • correspond to single- 
site case. 

Fig. 3: R(t) as a function of t in the cooling run for 
N = 256,(3 = 50,7 = 4,a = [-0.3,0.3] and j max ~ 1% 
of N 2 . The inset shows the fine structure of the peak. 

Fig. 4: A plot of D(Ra) versus Ra for the same pa- 
rameter values as in Fig. 3 for the multi-site case (o con- 
nected by dashed line) with a = 2.6. • with continuous 
line shows the single site case with a — 2.5. 
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